library(tidyverse)
library(parallel)
library(CoRC)
# helper to run tasks in parallel on all cores
mapInParallel <- function(data, fun, ..., .prep = {}) {
cl <- makeCluster(detectCores())
clusterCall(cl = cl, fun = eval, .prep, env = .GlobalEnv)
ret <- clusterApplyLB(cl = cl, x = imap(data, ~ set_names(list(.x), .y)), fun = lapply, as_mapper(fun), ...)
stopCluster(cl)
flatten(ret)
}This example loads the Kummer2000 - Oscillations in Calcium Signalling model. The model has 3 species which oscillate. These oscialltions can be visualized as a trajectory through a 3D space. The example does this once in a deterministic and once in a stochatic fashion.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000329")
#> # A copasi model reference:
#> Model name: "Kummer2000 - Oscillations in Calcium Signalling"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8
# Run 24 sec (2 Periods)
setTimeCourseSettings(24, intervals = 10000)
timeseries <- list(
deterministic = runTimeCourse()$result,
stochastic = runTimeCourse(method = list(method = "directMethod", use_random_seed = T, random_seed = 1))$result
)
# simplify species names so they are valid R syntax
timeseries <- map(timeseries, set_tidy_names, TRUE)
# Create the same plot for both timeseries
plots <- map(
timeseries,
plotly::plot_ly,
type = "scatter3d",
mode = "lines",
x = ~ G.alpha,
y = ~ activePLC,
z = ~ Calcium,
color = ~ Time
)
unloadModel()
plots$deterministicplots$stochasticThis implements an example from the Condor-COPASI paper. The example illustrates advantages of parallel processing.
# Run 1000 stochastic time series possibly in parallel
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
#> # A copasi model reference:
#> Model name: "Kummer calcium model"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8
# timeseries <- 1:1000 %>% map(~ runTimeCourse()$result)
timeseries <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
setTimeCourseSettings(method = list(method = "directMethod", use_random_seed = T))
}),
# iteration data (1000 random seeds),
1:1000,
# iteration code (formula ~)
~ runTimeCourse(method = list(random_seed = .x))$result
)
# Combine all results and reshape the data
plotdata <-
timeseries %>%
bind_rows() %>%
group_by(Time) %>%
# calculate mean and sd for all time points
summarise_all(funs(mean, sd)) %>%
# gather all values so the column "name" identifies "a_mean", "b_sd" etc.
gather("name", "value", -Time) %>%
# split up information on species (a,b,c) and type of value (mean, sd)
separate(name, c("species", "type"), "_") %>%
spread(type, value)
print(plotdata, n = 6)
#> # A tibble: 2,403 x 4
#> Time species mean sd
#> <dbl> <chr> <dbl> <dbl>
#> 1 0 a 8.00 0
#> 2 0 b 8.00 0
#> 3 0 c 8.00 0
#> 4 0.0500 a 7.06 0.249
#> 5 0.0500 b 8.12 0.117
#> 6 0.0500 c 5.60 0.442
#> # ... with 2,397 more rows
plot <-
ggplot(data = plotdata, aes(x = Time, y = mean, group = species, tt_sd = sd)) +
geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd, fill = species), alpha = 1 / 4) +
geom_line(aes(color = species)) +
guides(fill = "none") +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Concentration (", getQuantityUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("group", "x", "y", "tt_sd"))This implements an example from the Mendes2009 paper on COPASI use cases.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A copasi model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3
setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)
# Cartesian product of the input values
scan <- cross_df(
list(
cysteine = 0.3 * 10 ^ seq(0, 3, length.out = 6),
adomed = seq(0, 100, length.out = 51)
)
)
scan <-
scan %>%
mutate(
# Calculate steady state fluxes for every row
ss_fluxes = pmap(., function(cysteine, adomed) {
setSpecies("Cysteine", initial_concentration = cysteine)
setSpecies("adenosyl", initial_concentration = adomed)
ss <- runSteadyState()
stopifnot(ss$result == "found")
ss$reactions$flux
}),
# The second flux is CGS
CGS = map_dbl(ss_fluxes, 2),
# The third flux is TS
TS = map_dbl(ss_fluxes, 3)
)
plot <-
ggplot(data = scan, aes(x = adomed, group = cysteine)) +
geom_point(aes(y = CGS, color = "CGS")) +
geom_point(aes(y = TS, color = "TS")) +
geom_line(aes(y = CGS, color = "CGS")) +
geom_line(aes(y = TS, color = "TS")) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot)This implements an example from the Mendes2009 paper on COPASI use cases. It is in many ways similar to the previous example but is written to run parallelized.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A copasi model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3
# 10000 repeats of steady state task with random cysteine and adomed
scan <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)
}),
# iteration data (10000 random seeds),
1:10000,
# iteration code (formula ~)
~ {
set.seed(.x)
cysteine <- 0.3 * 10 ^ runif(1L, min = 0, max = 3)
adomed <- runif(1L, min = 0, max = 100)
setSpecies(
key = c("Cysteine", "adenosyl"),
initial_concentration = c(cysteine, adomed)
)
ss <- runSteadyState()
stopifnot(ss$result == "found")
list(
cysteine = cysteine,
adomed = adomed,
CGS = ss$reactions$flux[2],
TS = ss$reactions$flux[3]
)
}
)
# Combine all results and reshape the data
plotdata <-
scan %>%
bind_rows() %>%
gather("reaction", "flux", CGS, TS)
plot <-
ggplot(data = plotdata, aes(x = adomed, y = flux, group = reaction, tt_cys = cysteine)) +
geom_point(aes(color = reaction), alpha = 1 / 10, size = 3 / 4) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("tt_cys", "x", "y"))This implements an example from the Mendes2009 paper on COPASI use cases.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000010")
#> # A copasi model reference:
#> Model name: "Kholodenko2000 - Ultrasensitivity and negative feedback bring oscillations in MAPK cascade"
#> Number of compartments: 1
#> Number of species: 8
#> Number of reactions: 10
# get timeseries for the record
data_before <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
# read experimental data
data_experimental <-
read_tsv("data/MAPKdata.txt") %>%
rename(Time = time, "Mos-P" = "MAPKKK-P", "Erk2-P" = "MAPK-P") %>%
set_tidy_names(TRUE)
# define the experiments for copasi
fit_experiments <- defineExperiments(
data = data_experimental,
type = c("time", "dependent", "dependent"),
mapping = c(NA, "{[Mos-P]}", "{[Erk2-P]}"),
weight_method = "mean_square"
)
# define the parameters for copasi
fit_parameters <-
map(
parameter_strict(regex(c("V1$", "V2$", "V5$", "V6$", "V9$", "V10$"))),
~ {
val <- getParameters(.x)$value
defineParameterEstimationParameter(parameter(.x, "Value"), start_value = val, lower_bound = val * 0.1, upper_bound = val * 1.9)
}
)
result <-
runParameterEstimation(
parameters = fit_parameters,
experiments = fit_experiments,
method = "LevenbergMarquardt",
update_model = TRUE
)
# get timeseries for the record
data_after <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
plots <- list(
Erk2.P =
ggplot(mapping = aes(x = Time, y = Erk2.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Erk2-P (", getQuantityUnit(), ")"),
color = "Series"
),
Mos.P =
ggplot(mapping = aes(x = Time, y = Mos.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Mos-P (", getQuantityUnit(), ")"),
color = "Series"
)
)
unloadModel()
result$fitted_values
#> # A tibble: 2 x 5
#> fitted_value objective_value root_mean_square error_mean
#> <chr> <dbl> <dbl> <dbl>
#> 1 [Mos-P] 25.5 1.68 -0.283
#> 2 [Erk2-P] 10.1 1.06 0.330
#> # ... with 1 more variable: error_mean_std_deviation <dbl>
result$parameters
#> # A tibble: 6 x 8
#> parameter lower_bound start_value value upper_bound std_deviation
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 (MAPKKK activat… 0.250 2.46 2.46 4.75 0.139
#> 2 (MAPKKK inactiv… 0.0250 0.247 0.247 0.475 0.00828
#> 3 (dephosphorylat… 0.0750 0.882 0.882 1.42 0.261
#> 4 (dephosphorylat… 0.0750 1.42 1.42 1.42 0.709
#> 5 (dephosphorylat… 0.0500 0.728 0.728 0.950 0.0805
#> 6 (dephosphorylat… 0.0500 0.707 0.707 0.950 0.0962
#> # ... with 2 more variables: coeff_of_variation <dbl>, gradient <dbl>
result$protocol
#> [1] "Algorithm started at 2018-02-14 11:35:27.\nFor more information about this method see: http://www.copasi.org/tiki-index.php?page=OD.Levenberg.Marquardt\n\nIteration 379: Lambda reached max value. Terminating.\n\nAlgorithm reached the edge of the parameter domain 590 times.\n\nAlgorithm finished at 2018-02-14 11:35:29.\nTerminated after 380 of 2000 iterations.\n\n"
plotly::ggplotly(plots$Erk2.P)plotly::ggplotly(plots$Mos.P)